几种常见的判别分析分类方法在R中实现

判别分析（Discriminant analysis）用于基于一个或多个预测变量来预测对象属于某个给定类的概率，它对于连续或分类预测变量均适用。与logistic回归相比，判别分析在处理多分类问题时更稳健。

本篇汇总几种常见的判别分析类型的R操作过程，主要如下。

线性判别分析（Linear discriminant analysis，LDA），使用预测变量的线性组合预测给定的观测类别。LDA假设预测变量服从（多元）正态分布，并且类别间具有相同的方差（对于单变量分析）或相同的协方差矩阵（对于多变量分析）。

二次判别分析（Quadratic discriminant analysis，QDA），比LDA更灵活，假设预测变量服从（多元）正态分布，但允许各类别的协方差矩阵不等。

混合判别分析（Mixture discriminant analysis，MDA），假定每个类别都是子类别的高斯（正态）混合。

弹性判别分析（Flexible Discriminant Analysis ，FDA），使用预测变量的非线性组合，例如样条曲线（splines）。

正则化判别分析（Regularized discriminant anlysis，RDA），在预测变量数量大于训练集对象数量的情况下，正则化可改进协方差矩阵的估计。（注意其与冗余分析（Redundancy analysis，RDA）的缩写相同，不要混淆）

测试数据集

数据集概要

iris数据集，记录了150朵鸢尾花的花朵性状测量值。

这些鸢尾花来自三种物种，分别为setosa（n=50）、versicolor（n=50）和virginica（n=50）。

包含四种性状，分别为萼片长度（sepal length，cm）、萼片宽度（sepal width，cm）、花瓣长度（petal length，cm）和花瓣宽度（petal width，cm）。

#数据集，详情 ?iris
data(iris)
head(iris)

接下来期望从中找到合适的“变量组合”，作为区分不同鸢尾花的代表特征。 

*探索性分析

首先不妨查看一下各变量的数值分布，哪些花朵性状在不同物种之间具有明显的区别。

#变量分布
library(ggplot2)
 
ggplot(reshape2::melt(iris, id = 'Species'), aes(x = value)) +
geom_histogram(aes(fill = Species), color = 'gray') +
facet_wrap(~variable, ncol = 2, scales = 'free')

直方图清晰地表明，花瓣的长度和宽度似乎可以作为潜在特征，用以区分三种鸢尾花物种。

相比直接选取部分趋势明显的变量作为代表，通过降维技术所确定的变量组合特征通常是种更为实用的选择，例如主成分分析（PCA）。

#通过 PCA 降维数据，观测数据特征
pca <- princomp(iris[1:4])
plot(pca$scores[ ,1], pca$scores[ ,2], col = rep(c('red', 'green', 'blue'), summary(iris$Species)))

PCA显示，第一主成分可作为区分不同鸢尾花物种花朵性状的潜在特征。

接下来接入本篇的方法部分，结合已知的物种分类，使用带监督模式的判别分析实现数据降维，选择代表性的变量组合特征，以及实现分类。

线性判别分析（LDA）

LDA要求输入数据满足（多元）正态性，可通过QQ图评估。

#QQ-plot 检验多元正态性
qqplot(qchisq(ppoints(nrow(iris[1:4])), df = ncol(iris[1:4])), mahalanobis(iris[1:4], colMeans(iris[1:4]), cov(iris[1:4])))
abline(a = 0, b = 1)

QQ图显示了示例数据集满足多元正态性。

如果正态性假设被拒绝，可尝试转化数据的方式（如log转化，但要保证这种转化方式是合理的），获得正态分布的数据。

变量的标准化也是推荐的一步，对于消除变量间的量纲差异或者较大方差变量的权重时非常有效，可提高降维的精度。

#视情况选择标准化数据，如标准化为均值 0，标准差 1 的结构
iris[1:4] <- scale(iris[1:4], center = TRUE, scale = TRUE)

为了更好地展示LDA的分类器功能，将示例数据集分为两部分，一部分作为训练集用于LDA降维及分类器构建，另一部分作为测试集进一步评估LDA预测分类的功效。

#将数据集随机分为训练集（80%）和测试集（20%）
set.seed(123)
training <- sample(rownames(iris), nrow(iris)*0.8)
 
train.data <- subset(iris, rownames(iris) %in% training)
test.data <- subset(iris, ! rownames(iris) %in% training)

经过几个准备步骤后，执行LDA，这里以MASS包中的方法为例。

LDA算法首先查找使类别之间距离最大化的方向（即LDA轴），作为变量响应分类的最佳线性组合，并通过这种组合进一步预测分类。

library(MASS)
 
#拟合模型，详情 ?lda
model <- lda(Species~., data = train.data)
model

lda()确定各组数据的平均值并计算对象属于不同组的概率，将对象划分到概率最高的组中。

Prior probabilities of groups，各组的先验概率，即已知分组中所含对象数量占总数量的比例。例如在本示例中，随机抽取的训练集的setosa组中共含有40个对象（40个鸢尾花观测个体），占训练集所有对象（总计120个鸢尾花观测个体）的33.3%。

Group means，组均值，展示了每个分组中变量的平均值。

Coefficients of linear discriminants，线性判别系数，用于形成LDA决策规则的预测变量的线性组合。例如，LD1 = 0.828*Sepal.Length + 1.438*Sepal.Width - 2.179*Petal.Length - 2.656*Petal.Width。可在降维后根据线性判别系数所得表达式，预测训练集的分类。

Proportion of trace，可以将它理解为各LDA轴的“方差解释率”，以评估各LDA轴的重要程度，示例结果显示LDA第一主轴是非常重要的。

#作图观测对象在前两 LDA 轴中的投影
plot(model, col = rep(c('red', 'green', 'blue'), summary(train.data$Species)), dimen = 2)
 
#或者 ggplot2 作图
ggplot(cbind(train.data, predict(model)$x), aes(LD1, LD2, color = Species)) +
geom_point() +
stat_ellipse(level = 0.95, show.legend = FALSE)

结果显而易见了，LDA第一主轴承载了最能体现类间差异的特征。

有关ggplot2可视化的调整可见“排序图绘制”。

对比训练集中对象的既定分组属性和由LDA判别的分组属性的一致性，结果可表征LDA模型的拟合精度。

#对训练集预测分类
predictions <- predict(model, train.data)
 
#查看训练集对象的后验概率，即根据概率划分到高概率的类中
head(predictions$posterior)
#查看对训练集对象预测的分类
head(predictions$class)
#比较预测的分类和已知先验分类属性的差异，结果反映了准确度信息
mean(predictions$class == train.data$Species)

结果显示，98%以上的对象能够被分类到正确的类别中。

现在更改为测试集，进一步评估LDA分类器精度。

#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions$class == test.data$Species)

结果显示，约97%以的对象能够被分类到正确的类别中，LDA分类器的精度是相对可观的。

#后验概率也可通过热图展示
heatmap(predictions$posterior, Colv = NA, cexRow = 0.8, cexCol = 1)

实际应用中，还可以通过指定一个后验概率阈值，过滤模糊识别的结果，只保留较为可信的预测分类信息以提升优度等。

二次判别分析（QDA）

从某种意义上说，QDA由于不假定各类别的协方差相等，因此它比LDA灵活得多。这种优势主要体现在训练集数据量较大、或者观测类别较多时的情况，此时LDA要求的等协方差矩阵的假设经常被拒绝。

当训练集数据较少时，LDA往往比QDA更好。

MASS包也提供了QDA的方法。

#使用测试集拟合 QDA 模型，详情 ?qda
model <- qda(Species~., data = train.data)
model

qda()确定各组数据的平均值并计算对象属于不同组的概率，将对象划分到概率最高的组中。

Prior probabilities of groups，各组的先验概率，即已知分组中所含对象数量占总数量的比例。例如在本示例中，随机抽取的训练集的setosa组中共含有40个对象（40个鸢尾花观测个体），占训练集所有对象（总计120个鸢尾花观测个体）的33.3%。

Group means，组均值，展示了每个分组中变量的平均值。

随后，构建的分类器将对测试集对象预测分类。

#对训练集预测分类
predictions <- predict(model, train.data)
mean(predictions$class == train.data$Species)

结果显示，99%以上的训练集对象能够被分类到正确的类别中。

和上述LDA结果（使用相同的训练集数据，精度约~98%）相比，本示例也显示了QDA对训练集对象分类属性的识别有了提升，尽管甚微。

#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions$class == test.data$Species)

对于测试集，约97%以的对象能够被分类到正确的类别中，表明QDA分类器的精度也是可靠的。

混合判别分析（MDA）

无论LDA或QDA均要求每个类别的数据服从多元正态（高斯）分布。

对于MDA中存在的多个类别，均被假定为更多子类别的高斯混合，其中每个对象都有属于每个类别的概率。但MDA仍然假设各类别的协方差矩阵相等。

MDA可使用mda包实现。

library(mda)
 
#使用测试集拟合模型，详情 ?mda
model <- mda(Species~., data = train.data, subclasses = 3)
model
 
plot(model)

例如在上式中，定义每个既定类别均由3个子类别组成，通过这些子类别进行拟合分类。Training Misclassification Error代表了对训练集对象预测分类错误率，示例结果显示错误率很低。

散点图展示了降维后的数据特征，MDA第一主轴承载了最能体现类间差异的特征。

对于模型预测功能的实现，和上文方法类似。

#对训练集预测分类
predictions <- predict(model, train.data)
mean(predictions == train.data$Species) #结果等于 1-Training Misclassification Error

#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions == test.data$Species)

有些情况下，通过识别更小子类的MDA方法，可能会胜过全局识别模式LDA和QDA。

如下所示，存在3个主要的类别，每个类别又由3个子类别构成，各类别整体不满足正态（高斯）分布，但各子类别满足。由于子类别之间“不连续”，LDA和QDA的决策边界无法有效分类，但MDA可以通过正确识别子类别而实现良好的分类。

弹性判别分析（FDA）

FDA是LDA的灵活扩展，它使用样条曲线等预测变量的非线性组合。FDA可用于对每组内变量之间的多元非正态或非线性关系建模，从而实现更准确的分类。

FDA也可使用mda包实现。

#使用测试集拟合模型，详情 ?fda
model <- fda(Species~., data = train.data)
model

plot(model)

同样地，Training Misclassification Error代表了对训练集对象预测分类错误率，示例结果显示错误率很低。

散点图展示了降维后的数据特征，FDA第一主轴承载了最能体现类间差异的特征。

对于模型预测功能的实现，和上文方法类似。

#对训练集预测分类
predictions <- predict(model, train.data)
mean(predictions == train.data$Species) #结果等于 1-Training Misclassification Error

#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions == test.data$Species)

正则化判别分析（RDA）

RDA通过正则化组协方差矩阵建立分类规则，从而可以针对数据中的多重共线性提供更稳健的模型。这对于包含高度相关预测变量的大型多元数据集可能非常有用。

RDA可视为LDA和QDA之间的一种中间状态。RDA将QDA的独立协方差收缩为与LDA类似的等协方差。在预测变量的数量大于训练集对象数量的情况下，可以改善协方差矩阵的估计，潜在地提高了模型的准确性。 

RDA可使用klaR包实现。

library(klaR)

#使用测试集拟合模型，详情 ?rda
model <- rda(Species~., data = train.data)
model

Regularization parameters，gamma和lambda为两种正则化参数。

Prior probabilities of groups，各组的先验概率，即已知分组中所含对象数量占总数量的比例。例如在本示例中，随机抽取的训练集的setosa组中共含有40个对象（40个鸢尾花观测个体），占训练集所有对象（总计120个鸢尾花观测个体）的33.3%。

Misclassification rate，代表了对训练集对象预测分类错误率。

#对训练集预测分类
predictions <- predict(model, train.data)
mean(predictions$class == train.data$Species)

#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions$class == test.data$Species)

参考资料

http://www.sthda.com/english/articles/36-classification-methods-essentials/146-discriminant-analysis-essentials-in-r/

http://sebastianraschka.com/Articles/2014_python_lda.html

友情链接